Woodbury Optimizations by aschneuw · Pull Request #4 · r-gam/gamm4

aschneuw · 2026-04-06T16:39:35Z

Improves Woodbury Identity optimizations
properly takes into account pivoting for all Cholesky Decomposition based operations
updates scikit-sparse to the latest version 0.5.0 (had breaking API changes)
Introduces Ztilde_t <- phi.rphi %*% Zt as a further optimization

aschneuw · 2026-04-08T09:06:32Z

    if (length(ind)) { ## extract columns corresponding to non-smooth r.e.s 
      Zt <- getME(ret$mer,"Zt")[ind,] ## extracting random effects model matrix
-      root.phi <- getME(ret$mer,"Lambdat")[,ind] ## and corresponding sqrt of cov matrix (phi) (was ind,ind - seems wrong)
+      root.phi <- getME(ret$mer,"Lambdat")[ind,ind] ## and corresponding sqrt of cov matrix (phi) (was ind,ind - seems wrong)


this was exactly my thinking :)

fabian-s

Review Summary

Thanks for the work on this, Arno — the getVb() refactoring and the Woodbury identity path are well-structured and the core linear algebra is correct. Numerical equivalence testing (7 scenarios: Gaussian/binomial/Poisson with various smooth types, random intercepts/slopes) confirms machine-precision agreement between the R-only path and master, and between the Woodbury and direct paths.

However, there are several issues that need to be addressed before this can be merged. See inline comments for details.

Must Fix

python_cholmod hardcoded to TRUE — will crash any user without Python/scikit-sparse
Permutation direction in phi.rphi — likely uses perm where inv_perm is needed
root.phi potentially undefined — when no non-smooth random effects exist

Should Fix

Comment typo introduced ("--get need")
Stale/misleading comment on root.phi line
scikit-sparse==0.5.0 exact pin is fragile
Dead V computation
Dead if (TRUE) {} else {} branches
Indentation inconsistency

fabian-s · 2026-04-08T10:22:48Z

    if (length(ind)) { ## extract columns corresponding to non-smooth r.e.s 
      Zt <- getME(ret$mer,"Zt")[ind,] ## extracting random effects model matrix
-      root.phi <- getME(ret$mer,"Lambdat")[,ind] ## and corresponding sqrt of cov matrix (phi) (was ind,ind - seems wrong)
+      root.phi <- getME(ret$mer,"Lambdat")[ind,ind] ## and corresponding sqrt of cov matrix (phi) (was ind,ind - seems wrong)


Stale comment: (was ind,ind - seems wrong) — but the code IS using [ind,ind], same as master. The parenthetical should be removed to avoid confusion.

fabian-s · 2026-04-08T10:26:36Z

Sorry for the chaos - briefly lost control of my coding agent ... ;)
I fully endorse the comments it made on my behalf above.

fabian-s · 2026-04-13T08:54:34Z

I can’t push to aschneuw/r-gam-gamm4 from here (permission error on the fork), so I’m putting the proposed fix set here directly, see also https://github.com/r-gam/gamm4/tree/woodbury-fixes which you could simply merge into your PR.

I rechecked this locally with sksparse enabled via reticulate, and I think the remaining changes needed on this branch are:

Restore the Python Woodbury fallback when factoring phi fails.
Fix the direct python_cholmod path: the current manual permutation handling is not correct for scikit-sparse 0.5.0, and the current V_py = r_to_py(V) also hits a read-only buffer failure.
Reset the accidental executable bits on the touched package files.
Keep the small cleanup items (--get need typo, dead V <- Diagonal(...) line).

The direct-path fix is to use cho_factor(..., lower=TRUE) plus solve(system = "A"), because in sksparse 0.5.0 that is the solve mode that applies CHOLMOD’s permutation correctly.

Suggested code changes in R/gamm4.r:

if (woodbury && nrow(Zt)) {
  if (python_cholmod) {
    phi_py = reticulate::r_to_py(phi)$copy()
    cho_phi_py <- try(
      cholmod$cho_factor(phi_py, supernodal_mode="auto", lower=TRUE),
      silent = TRUE
    )
    if (inherits(cho_phi_py, "try-error")) {
      woodbury <- FALSE
    } else {
      perm <- as.integer(reticulate::py_to_r(cho_phi_py$get_perm()) + 1)
      inv_perm <- perm
      inv_perm[perm] <- 1:length(perm)
      phi.rphi <- reticulate::py_to_r(scipy_sparse$csc_matrix(cho_phi_py$R))[,inv_perm]
    }
  } else {
    phi.inv <- try(chol2inv(chol(phi)), silent = TRUE)
    if (inherits(phi.inv,"try-error")) woodbury <- FALSE
  }
}

if (python_cholmod) {
  # Only system="A" applies CHOLMOD's internal permutation correctly.
  V_py = reticulate::r_to_py(V)$copy()
  R <- cholmod$cho_factor(V_py, supernodal_mode="auto", lower=TRUE)
  XFP_py <- reticulate::r_to_py(Xfp)$copy()
  XF_py <- reticulate::r_to_py(Xf)$copy()
  XVXS <- t(Xfp) %*% reticulate::py_to_r(scipy_sparse$csc_matrix(R$solve(XFP_py, system="A"))) + Sp^2/scale
  XVX <- t(Xf) %*% reticulate::py_to_r(scipy_sparse$csc_matrix(R$solve(XF_py, system="A")))
} else {
  R <- mgcv::mchol(V);piv <- attr(R,"pivot")
  WX <- as(solve(t(R),Xfp[piv,]),"matrix")
  XVX <- as(solve(t(R),Xf[piv,]),"matrix")
  XVX <- crossprod(XVX)
  XVXS <- crossprod(WX)+Sp^2/scale
}

gam.terms <- attr(gmf,"terms") # terms object for `gam' part of fit -- need this for prediction to work properly

and remove the dead line:

V <- Diagonal(n=length(v),x=v)

I also recommend resetting the file modes on:

.gitignore
ChangeLog
DESCRIPTION
MD5
NAMESPACE
R/gamm4.r
man/gamm4.Rd

Locally I validated the patched version with:

python_cholmod=TRUE vs FALSE model-level smoke fits (Gaussian smooth-only, Gaussian + random intercept, Poisson + random intercept): all matched to numerical precision
randomized direct-path getVb(..., woodbury=FALSE) comparisons against the native-R branch: no mismatches after the direct-path fix above
a singular / semidefinite phi case: the Woodbury Python branch now falls back cleanly instead of hard-failing during factorization

From my side, once the above is applied, this looks ready to merge.

… and further optimizations

aschneuw · 2026-04-13T20:10:52Z

I can’t push to aschneuw/r-gam-gamm4 from here (permission error on the fork), so I’m putting the proposed fix set here directly, see also https://github.com/r-gam/gamm4/tree/woodbury-fixes which you could simply merge into your PR.

I rechecked this locally with sksparse enabled via reticulate, and I think the remaining changes needed on this branch are:

Restore the Python Woodbury fallback when factoring phi fails.

Fix the direct python_cholmod path: the current manual permutation handling is not correct for scikit-sparse 0.5.0, and the current V_py = r_to_py(V) also hits a read-only buffer failure.

Reset the accidental executable bits on the touched package files.

Keep the small cleanup items (--get need typo, dead V <- Diagonal(...) line).

The direct-path fix is to use cho_factor(..., lower=TRUE) plus solve(system = "A"), because in sksparse 0.5.0 that is the solve mode that applies CHOLMOD’s permutation correctly.

Suggested code changes in R/gamm4.r:
if (woodbury && nrow(Zt)) {
  if (python_cholmod) {
    phi_py = reticulate::r_to_py(phi)$copy()
    cho_phi_py <- try(
      cholmod$cho_factor(phi_py, supernodal_mode="auto", lower=TRUE),
      silent = TRUE
    )
    if (inherits(cho_phi_py, "try-error")) {
      woodbury <- FALSE
    } else {
      perm <- as.integer(reticulate::py_to_r(cho_phi_py$get_perm()) + 1)
      inv_perm <- perm
      inv_perm[perm] <- 1:length(perm)
      phi.rphi <- reticulate::py_to_r(scipy_sparse$csc_matrix(cho_phi_py$R))[,inv_perm]
    }
  } else {
    phi.inv <- try(chol2inv(chol(phi)), silent = TRUE)
    if (inherits(phi.inv,"try-error")) woodbury <- FALSE
  }
}
if (python_cholmod) {
  # Only system="A" applies CHOLMOD's internal permutation correctly.
  V_py = reticulate::r_to_py(V)$copy()
  R <- cholmod$cho_factor(V_py, supernodal_mode="auto", lower=TRUE)
  XFP_py <- reticulate::r_to_py(Xfp)$copy()
  XF_py <- reticulate::r_to_py(Xf)$copy()
  XVXS <- t(Xfp) %*% reticulate::py_to_r(scipy_sparse$csc_matrix(R$solve(XFP_py, system="A"))) + Sp^2/scale
  XVX <- t(Xf) %*% reticulate::py_to_r(scipy_sparse$csc_matrix(R$solve(XF_py, system="A")))
} else {
  R <- mgcv::mchol(V);piv <- attr(R,"pivot")
  WX <- as(solve(t(R),Xfp[piv,]),"matrix")
  XVX <- as(solve(t(R),Xf[piv,]),"matrix")
  XVX <- crossprod(XVX)
  XVXS <- crossprod(WX)+Sp^2/scale
}
gam.terms <- attr(gmf,"terms") # terms object for `gam' part of fit -- need this for prediction to work properly
and remove the dead line:
V <- Diagonal(n=length(v),x=v)
I also recommend resetting the file modes on:

.gitignore

ChangeLog

DESCRIPTION

MD5

NAMESPACE

R/gamm4.r

man/gamm4.Rd

Locally I validated the patched version with:

python_cholmod=TRUE vs FALSE model-level smoke fits (Gaussian smooth-only, Gaussian + random intercept, Poisson + random intercept): all matched to numerical precision

randomized direct-path getVb(..., woodbury=FALSE) comparisons against the native-R branch: no mismatches after the direct-path fix above

a singular / semidefinite phi case: the Woodbury Python branch now falls back cleanly instead of hard-failing during factorization

From my side, once the above is applied, this looks ready to merge.

Applied your suggested changes. Removed "# Only system="A" applies CHOLMOD's internal permutation correctly." Since this was misleading to me. I think the previous approach would have been doable as well with the appropriate permutations applied. But your approach works as well and is better.

fabian-s

please resolve my earlier comment on line gamm4.r l. 587
since the python_cholmod = TRUE branch needs py_require and that was only added in reticulate 1.41, please add the version requirement to DESCRIPTION

fabian-s · 2026-05-04T09:25:56Z

@simonnwood this is now ready to merge from my perspective.
please let us know if you require additional changes.

aschneuw commented Apr 8, 2026

View reviewed changes

fabian-s reviewed Apr 8, 2026

View reviewed changes

Comment thread R/gamm4.r Outdated

fabian-s reviewed Apr 8, 2026

View reviewed changes

Comment thread R/gamm4.r Outdated

fabian-s reviewed Apr 8, 2026

View reviewed changes

Comment thread R/gamm4.r

fabian-s reviewed Apr 8, 2026

View reviewed changes

Comment thread R/gamm4.r Outdated

fabian-s requested changes Apr 8, 2026

View reviewed changes

aschneuw requested a review from fabian-s April 10, 2026 20:55

aschneuw force-pushed the woodbury-fixes branch from a89eb80 to b1f6854 Compare April 12, 2026 17:30

aschneuw force-pushed the woodbury-fixes branch 2 times, most recently from 86d3303 to 1754ea1 Compare April 13, 2026 19:47

Update scikit-sparse to 0.5.0, pivoting for Python cholmod operations…

9574240

… and further optimizations

aschneuw force-pushed the woodbury-fixes branch from 1754ea1 to 9574240 Compare April 13, 2026 19:53

fabian-s requested changes Apr 28, 2026

View reviewed changes

aschneuw added 2 commits April 28, 2026 13:51

Min reticulate 1.41.0

85b06c7

Fix undefined root.phi

ee08496

aschneuw requested a review from fabian-s April 30, 2026 13:49

fabian-s approved these changes May 4, 2026

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Woodbury Optimizations#4

Woodbury Optimizations#4
aschneuw wants to merge 3 commits into
r-gam:masterfrom
aschneuw:woodbury-fixes

aschneuw commented Apr 6, 2026

Uh oh!

aschneuw Apr 8, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

fabian-s left a comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

fabian-s Apr 8, 2026

Uh oh!

Uh oh!

fabian-s commented Apr 8, 2026

Uh oh!

fabian-s commented Apr 13, 2026 •

edited

Loading

Uh oh!

aschneuw commented Apr 13, 2026

Uh oh!

fabian-s left a comment

Uh oh!

fabian-s commented May 4, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

aschneuw commented Apr 6, 2026

Uh oh!

aschneuw Apr 8, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

fabian-s left a comment

Choose a reason for hiding this comment

Review Summary

Must Fix

Should Fix

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

fabian-s Apr 8, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

fabian-s commented Apr 8, 2026

Uh oh!

fabian-s commented Apr 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

aschneuw commented Apr 13, 2026

Uh oh!

fabian-s left a comment

Choose a reason for hiding this comment

Uh oh!

fabian-s commented May 4, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

fabian-s commented Apr 13, 2026 •

edited

Loading